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ABSTRACT 


Time-frequency signal representations combined with subspace identification methods were used to analyze 
aeroelastic flight data from the F/A- 1 8 Systems Research Aircraft (SRA) and aeroservoelastic data from the F/A- 1 8 
High Alpha Research Vehicle (HARV). The F/A- 18 SRA data were produced from a wingtip excitation system that 
generated linear frequency chirps and logarithmic sweeps. HARV data were acquired from digital Schroeder-phased 
and sine pulse excitation signals to actuator commands. Nondilated continuous Morlet wavelets implemented as a 
filter bank were chosen for the time-frequency analysis to eliminate phase distortion as it occurs with sliding window 
discrete Fourier transform techniques. Wavelet coefficients were filtered to reduce effects of noise and nonlinear 
distortions identically in all inputs and outputs. Cleaned reconstructed time domain signals were used to compute 
improved transfer functions. Time and frequency domain subspace identification methods were applied to enhanced 
reconstructed time domain data and improved transfer functions, respectively. Time domain subspace performed 
poorly, even with the enhanced data, compared with frequency domain techniques. A frequency domain subspace 
method is shown to produce better results with the data processed using the Morlet time-frequency technique. 

NOMENCLATURE 


a 

wavelet scaling parameter 

[AJi,C,D] 

system state-space matrices 

AOA 

angle of attack 

ASE 

aeroservoelastic 

CWT 

continuous wavelet coefficient 

da 

finite interval of a 

DFT 

discrete Fourier transform 

e 

exponential function 

f 

center frequency of wavelet filter 

SF 

Fourier transform 

FIR 

finite impulse response 

G 

transfer function 

h 

wavelet basis function 

H 

Toeplitz matrix of Markov parameters 

//(CO) 

Fourier transform of wavelet basis function 

H f 

wavelet filter centered at/ 

HARV 

High Alpha Research Vehicle 

I 

identify matrix 

j 

imaginary J~\ 

K 

Kalman gain matrix 

M 

number of frequency data 

MIMO 

multi-input-multi-output 

n 

number of system states 

N 

number of time data 

OBES 

on-board excitation system 

P 

number system outputs 



[P,Q,S] 

matrices 

SISO 

single-input-single-output 

SRA 

Systems Research Aircraft 

t 

time 

T 

total time of wavelet analysis data 

u 

system input vector 

U 

block Hankel input matrix 

V 

block Hankel noise matrix 

X 

system state vector 

X 

state matrix 

x(t) 

analysis data in time domain 

X«o) 

analysis data in frequency domain 

y 

system output vector 

Y 

block Hankel output matrix 

z 

discrete frequency 

r 

extended observability matrix 

X 

wavelet time shift parameter 

0) 

radian frequency 

Superscripts 

a 

antisymmetric 

S 

symmetric 

Subscripts 

i,k 

integer indices 

o 

noise subspace 

s 

signal subspace 

Symbols 

_L 

orthogonal complement 

t 

pseudoinverse 

A 

estimate 

U 

omission of first rows of matrix U 

U 

omission of last rows of matrix U 


INTRODUCTION 

Envelope expansion of new or modified aircraft often necessitates structural stability testing to verify safety 
margins and to prevent catastrophic flutter. Inflight testing allows determination of aeroelastic or aeroservoelastic 
effects as a function of flight parameters. Flight data are acquired so stability estimation and system identification 
can be compared with analytic predictions. Any anomalies are regarded with care to guarantee maximum flight 
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safety. Because atmospheric turbulence is generally insufficient to determine modal characteristics, excitation 
systems are often essential to resolve stability trends from noisy measurements. 

Flight flutter data processing attempts to analyze data characteristics in terms of timing, shape, amplitude, 
frequency, and duration of events in the data. These analyses have traditionally been addressed using classical 
Fourier techniques. However, Fourier methods are suspect due to the inherently transient nature of inflight 
aeroelastic dynamics. The infinite and at least locally-periodic waveform assumptions in Fourier analysis cannot 
adequately describe the intermittency, modulation (amplitude, phase, or frequency), nonperiodicity, non- 
stationariness, time- variance, or nonlinearity in the data. 

This paper describes a novel filtering procedure which removes distortions from aeroelastic or aeroservoelastic 
(ASE) flight data. The filtering procedure builds upon and adapts the methods of wavelet transform analysis to the 
particular test data generated from structural excitation systems. 

Wavelets are versatile harmonic analysis tools which combine time and frequency representations into localized 
waveforms. Given a segment of experimental data, the wavelet transform convolves a selected series of local 
waveforms with the data to identify correlated features, or patterns, in the signal. The resulting set of wavelet 
coefficients can be interpreted as multidimensional correlation coefficients. Features of shape, size, and location are 
naturally characterized by these waveforms and related coefficients. 

The salient features of the original signal may be reconstructed by exploiting the redundancy of the wavelets in 
the continuous wavelet transform. Time-frequency components can be removed from these data by masking off the 
unwanted components (setting the corresponding wavelet coefficients to zero). This method is used to filter 
unwanted distortions and extract desired features from the input (excitation) and output (structural response) data. 
This feature extraction method offers advantages to traditional band-pass filtering or thresholding techniques in that 
it removes unwanted features while leaving desirable components intact. 

By removing aspects of system inputs and outputs which are detrimental to linear identification methods, 
improvement of system identification is realized. Subspace identification methods have recently shown promise as 
being efficient, generally applicable, and robust. Experience with these methods will be discussed and comparisons 
made between transfer functions estimated with and without feature processing. 

AIRCRAFT DESCRIPTIONS AND TEST PROCEDURES 

Descriptions of the F/A-18 SRA and HARV follow. Flight flutter testing using a wingtip excitation system on 
the SRA is described in detail. High angle-of-attack aeroservoelastic research with the HARV is described using 
commands generated from the control system. 

Systems Research Aircraft 


The F/A-18 Systems Research Aircraft (SRA) is being flight tested at NASA Dryden Flight Research Center, 
Edwards, California, for primarily flight systems experiments. These tests include optical sensing, new actuation 
concepts, smart structures, and advanced airdata and flight control systems. A major left wing structural modification 
was done, however, on the SRA to allow testing of several hydraulic and electromechanical advanced aileron 
actuator concepts. Because the test actuators may be larger than the standard one, a fitting called a hinge-half 
(supporting the aileron hinge, the actuator, and a fairing) had to be replaced by larger and heavier items. 
Approximately 35 lb were added to the wing. 

These SRA data were derived from flight flutter testing using a wingtip aerodynamic structural excitation 
system 1 at subsonic, transonic, and supersonic flight conditions. 2 Inputs were linear sinusoidal frequency sweeps, or 
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chirps, and logarithmic sweeps of the excitation system. Each wingtip exciter had a small fixed aerodynamic vane 
forward of a rotating slotted hollow cylinder (fig. 1). When the cylinder rotated, the aerodynamic pressure 
distribution on the vane generated a force. This force changed at twice the cylinder rotation frequency. The rotation 
direction of the cylinder could generate two force levels, low-force levels (at 25 percent open) or high-force levels 
(at 75 percent open). Exciters could be placed in various locations. Three possibilities were mid-wingtip, forward 
wingtip at the wingtip leading edge, or aft wingtip near the wingtip trailing edge (as pictured in figure 1). 

Table 1 summarizes the matrix of test points for the SRA. Wingtip exciters were mounted on the SRA wingtip 
launcher rails in forward and aft positions, independently. They were operated in both low- and high-force modes 
for 15-, 30-, or 60-second forward and reverse frequency sweeps. A total of 260 test points were flown with exciters 
in one or more of the optional configurations at specified flight conditions. Linear and logarithmic sinusoidal sweeps 
of the exciter up to 40 Hz were used to excite the primary modes of interest in the 5 to 30 Hz range (table 2). Some 
sweeps were actually multiples of two or four shorter duration sweeps with no interrupt such as double (sweep up, 
then down) and quadruple (up-down-up-down) sweeps. Generally sweeps were performed for symmetric and 
antisymmetric excitation in each configuration. 



Figure 1 . Wingtip exciter on F/A-18 SRA. 


Table 1 . SRA aeroelastic flight test matrix. 


Mach number 
Altitude, ft 
Exciter position 
Force 
Sweep 

Duration, sec 
Range, Hz 
Multiple 


0.54, 0.65, 0.70, 0.80, 0.85, 0.90, 0.95, 1.05, 1.2, 1.4, 1.6 
10k, 30k, and 40k 

Both forward, both aft, and left-aft/right-forward 

Low and high 

Linear and logarithmic 

15, 30, and 60 

3-12, 25-35, 3-35, 35-3, and 3^10 

Single, double, and quadruple contiguous sweeps per maneuver 




Table 2. F/A-18 SRA with wingtip exciter: Calculated elastic modal frequencies. 


Symmetric Mode 

Hz 

Antisymmetric Mode 

Hz 

Wing first bending 

5.59 

Fuselage first bending 

8.15 

Fuselage first bending 

9.30 

Wing first bending 

8.84 

Stabilizer first bending 

13.21 

Stabilizer first bending 

12.88 

Wing first torsion 

13.98 

Wing first torsion 

14.85 

Fin first bending 

16.83 

Fin first bending 

15.61 

Wing second bending 

16.95 

Wing second bending 

16.79 

Wing outboard torsion 

17.22 

Fuselage second bending 

18.62 

Fuselage second bending 

19.81 

Trailing-edge flap rotation 

23.47 

Trailing-edge flap rotation 

23.70 

Fuselage torsion 

24.19 

Stabilizer fore and aft 

28.31 

Launcher rail lateral 

24.35 

Wing second torsion 

29.88 

Stabilizer fore and aft 

28.58 

Aileron rotation 

33.44 

Wing second torsion 

29.93 

Aileron torsion 

38.60 

Aft fuselage torsion 

37.80 

Stab second bending, wing third bending 

43.17 

Wing pitch 

39.18 


High Alpha Research Vehicle (HARV) 

Another F/A-18 aircraft (fig. 2) was modified at the NASA Dryden Flight Research Center to perform flight 
research at high angle of attack (AOA) using thrust vectoring and incorporating control law concepts for agility and 
performance enhancement, as well as providing data for correlation with computational fluid dynamics solutions. 
This vehicle is referred to as the High Alpha Research Vehicle (HARV) 3 throughout this report. As opposed to the 
SRA, the HARV was structurally very different from a standard F/A-18 aircraft, as can be seen by comparing the 
HARV modal data in table 3 with the SRA data in table 2. 

Structural modifications included addition of Inconel® vanes in each engine exhaust for thrust vectoring. 
Corresponding ballast was added in the forward fuselage to maintain the center-of-gravity location. A research 
flight control system was incorporated for feedback control of aerodynamic surfaces and the vanes. An inertial 
navigation system was installed for AOA and angle of sideslip rate feedbacks, wingtip AOA vanes, and pressure 
probes (for airdata research purposes). Additional instrumentation was added for loads, vane temperatures, and 
structural dynamics. 4 

An important element of the HARV flight system was the onboard excitation system (OBES). This system was 
implemented to add programmed digital signals to the control system actuator commands for structural excitation. 
Inputs from 5 to 25 Hz were added to the stabilator, aileron, rudder, and pitch- and yaw-vectored thrust commands. 
Data were generated with OBES commands at 5° to 70° AOA 4 at 1 g. Accelerometers were located in the aircraft 
nose, vertical and horizontal tails, and wingtips. 


®Inconel is a registered trademark of Huntington Alloy Products Div., International Nickel Co., Huntington, West Virginia. 
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Figure 2. F/A-18 HARV. 


Table 3. F/A-18 HARV calculated elastic modal frequencies. 


Symmetric Mode 

Hz 

Antisymmetric Mode 

Hz 

Wing first bending 

5.79 

Fuselage first bending 

7.13 

Fuselage first bending 

7.71 

Wing first bending 

8.75 

Wing first torsion 

11.68 

Wing first torsion 

12.03 

Stabilizer first bending 

13.73 

Stabilizer first bending 

13.63 

Wing fore-aft 

18.53 

Wing fore-aft 

15.18 

Fin first bending 

15.92 

Fin first bending 

15.71 

Wing second bending 

17.55 

Fuselage first torsion 

19.13 

Fuselage second bending 

15.87 

Fuselage second bending 

21.42 

Exhaust vane rotation 

22.10 

Exhaust vane rotation 

22.10 

Inboard flap rotation 

23.60 

Inboard flap rotation 

23.18 

Wing outboard torsion 

27.50 

Fore-fuselage torsion 

24.17 


HARV data were generated from digital sine j pulses and Schroeder-phased 56 inputs into control surface 

commands at generally high AOA. Sine pulses and Schroeder-phased signals have a theoretically flat frequency 
response across a defined frequency band and are optimal for structural excitation because all modes are excited 
equally. Furthermore, Schroeder-phased signals have a minimal peak crest factor (difference between upper and 
lower bound of the power spectrum) amongst all signals with the same total power and are therefore preferred. 
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Schroeder-phased signals are multifrequency signals composed of a large number of harmonics equally spaced 
in frequency. Each harmonic is specified by a phase shift so that when they are summed the waveform has a low peak 
factor and fits a given power spectrum. It is represented by 

x(t) = A^cos(wr + ^-) (1) 

k 

where N is the number of data, and A is a specified amplitude to acquire maximum power under the saturation limits. 

TIME-FREQUENCY ANALYSIS METHOD 

This section illustrates the limitations of the Fourier transform as a tool for analysis of transient or short time 
period phenomena. The wavelet transform is offered as an alternative analysis tool, and the Morlet wavelet is 
selected as the wavelet basis function. Criteria for selecting the frequency width of the transform are discussed. 
Finally, a method for reconstructing the desirable features of the original signal from the wavelet coefficients and 
Morlet wavelet basis function is offered. 


Wavelet Transform 

Analysis of the frequency content of a signal using the Fourier transform results in a spectral representation that 
is a function of frequency. This transform is not localized in time. For a time-varying signal x(t), however, a 
transform of the form 


T(x(t)) = X(t, CO) (2) 

is required to locate instantaneous frequencies. Such a transform is limited in resolution by the Heisenberg 
uncertainty principle 


AcoAr>2 (3) 

Resolution problems resulting from the uncertainty principle can be minimized by extending the idea of the 
Fourier transform to a new time-frequency decomposition called the wavelet decomposition. 1 This decomposition 
performs an orthonormal projection of the signal onto a set of basis functions that are adapted to the required 
frequency resolution. The Continuous Wavelet Transform ( CWT) is defined as 

CWT(X,a) = f 00 x(.t)h *(— 1 dt (4) 

where h a T (r) = h (j— j is the wavelet basis function, X is the local time, and a is a scale or dilation parameter set 

to match the level of resolution desired. Wavelet transforms can be looked upon as the filtering of a signal through 
a bank of filters h a T (r). These filters are band-pass with identical shape but with varying frequency width da centered 
around the frequencies of interest. Note that the h z (t) are not orthogonal because they are redundant (defined for 
continuous a and x). With the CWT the time resolution is arbitrarily good at high frequencies, and the frequency 
resolution becomes arbitrarily good at low frequencies (within the limits of the uncertainty principle) because da/ a 
is constant. Time-frequency analyses attempt to minimize the effects of these resolution problems by using a priori 
knowledge of the signal properties or adapting the resolution to the signal. 8 ’ 9 
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Selection of Wavelet Basis 


Aeroelastic and aeroservoelastic excitation data are essentially short-time sinusoidal, so the selection of a 
wavelet basis function should encompass this characteristic. Therefore, the Morlet wavelet 10 


h{t) = 




(5) 


was chosen as the basis function because of its clear interpretation in the frequency domain (Gaussian window) and 
time domain (locally periodic waveform) for the analysis of vibration data. Figures 3(a) and 3(b) show real and 
imaginary components of the Morlet wavelet. Figure 3(c) shows a Fourier representation of the real part. 




Time, sec 


970002 


(b) Imaginary part (middle). 
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(c) Fourier transform of real part (bottom). 
Figure 3. Morlet continuous wavelet transform. 



For analysis of the SRA data composed of linear frequency sweeps where the rate of change is constant over the 
whole time interval, a constant resolution is required along the entire time-frequency domain. Hence, the chosen 
time-frequency transform uses nondilated ( a = constant) Morlet wavelets as functions onto which the data are 
projected. For a selected frequency and frequency resolution, the wavelet transform convolves the corresponding 
Morlet filter with the signal to produce the wavelet coefficients. 

Resolution Criteria 


A tradeoff between time and frequency resolution is used to select the appropriate frequency resolution of the 
filter. In the case of the sweeps considered here, the frequency varies from to, to tt^ in T seconds. Assuming a desired 
resolution of Ago in frequency, the time resolution Ar is actually given by the time localization of this frequency by 

Ar = — — — Ago (6) 

go 2 - W] 


From the Heisenberg uncertainty principle (eq. (3)), we have 


Aco 


2 



and At 



(7) 


For a typical linear sweep from 0 to 30 Hz in 30 sec, these relations give A t = Jl sec, and Aco = J2. Hz. These 
values will be used to determine the frequency width of the Morlet filters. 


Time-Frequency Filtering and Signal Reconstruction 


In the frequency domain, the wavelet coefficient corresponding to a Morlet filter with center frequency /is 

CWT(co, /) = 3T(CWT(t, /)) = X(co) • H f (t o) (8) 


where the Fourier transform of the Morlet wavelet is 


H f ( to) s //(go, /) (9) 

The complex Morlet basis function from equation (5) presents problems because it should have no group delay 
(real Fourier representation) and no phase distortion. For example, when reconstructing the original time domain 
signal from a bank of complex Morlet filters, phase distortion results. Therefore, only the real part of the Morlet 
wavelet is used to ensure that the transform HJ<0) is real, and signal reconstruction is in-phase with the original 
signal. The Fourier transform of the real part 


H f ( to) = e 


-(to-/) 


- e 


2 -2 
-03 -/ 


( 10 ) 


is shown in figure 3(c). Convolution and inverse Fourier transform operations produce the desired time- 
frequency coefficients 

CWT(t,f) = ^■ 1 (X(to)// / (GD)) (11) 


where X(co) and HA go) are the Fourier transforms of signal and Morlet wavelet, respectively. This real Morlet filter 
bank provides a finite impulse response (FIR) linear phase filter for each center frequency / and maintains phase 
consistency between the original signal and some type of reconstruction. 
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Assuming that n wavelet basis functions are to be used to perform the signal reconstmction, equation (8) is used 
to construct the FIR filter bank. For a given frequency co 0 , the wavelet coefficients are related to the original signal 
spectrum by 


CWT( co 0> /,) 


*/,(««>> 


= X(to 0 ) 


CWT((O 0 ,f n )_ 


H f (®q) 

•> n 


This equation is solved for X(co Q ) using the pseudoinverse at every frequency co 0 


X(w 0 ) 


CWT(<Hq, / j) 

’"/,(» o) 

CWT((a 0 , /„) 

H f (0) 0 ) 

n 


(13) 


and is inverse transformed to get the time domain reconstruction x(t) . 

This procedure will be used subsequently to filter unwanted distortions and extract desired features from the 
input excitation and output structural response data from the SRA and HARV aircraft. These feature-filtered data 
will then be used to compute cleaned signals for transfer functions and subspace identification algorithms. Figure 4 
shows an outline of the entire procedure. 

SIGNAL PROCESSING IN THE TIME-FREQUENCY DOMAIN 

Transfer functions are routinely used in aeroelastic and aeroservoelastic analyses to acquire state-space 
representations of the system modal dynamics and to predict flutter boundaries with damping trends or other 
advanced techniques. 11 - 12 Traditional Fourier transform methods often disguise important features in the data from 



Figure 4. Time-frequency Morlet feature-filtered system identification (ID) procedure. 
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averaging and windowing. To circumvent these problems, a recipe using the wavelet-based feature extraction filter 
is employed to estimate cleaner transfer functions based on time-frequency localization. Figure 4 shows a schematic 
of the entire procedure. 

• Feature filter by masking the dominant harmonic band in the inputs and performing the same procedure to 
the outputs. The resulting input-output pairs are from a strictly linear system. 

• Localize the relevant frequency content by filtering each input-output signal with a bank of narrow 
band-pass FIR filters centered around each co ( . 

• Localize the relevant time information by windowing each band-pass-filtered signal around the time point 
where a filtered input has maximum amplitude. The relevant information at frequency co ( is assumed to be 
localized in time. Energy outside the chosen time window at identical frequencies is considered as noise. 

• Apply equation (13) for every windowed input-output pair at each C0q. 

• Reconstruct time signals x(t) for time domain identification or derive transfer functions using the spectral 
input-output pairs of X(to) . 

Figure 5 shows an example of what the Morlet filter can do beyond standard Fourier techniques. The top four 
plots are derived from standard Fourier techniques, from one exciter (left) to accelerometer output and other exciter 
(right) to the same output. The bottom four plots are the corresponding Morlet feature-filtered results. Improvement 
is obvious with the Morlet-processed data in identifying modal peaks from the noise and establishing well-defined 
phase response. 

A comparison is made between the Morlet wavelet filter and signal analysis with a discrete Fourier transform 
(DFT) using a sliding window in figure 6. The signal is a linear sweep measured from the SRA exciter vane as the 
aeroelastic system input. The sliding DFT is computed for the same frequency resolution as the Morlet filter bank. 
A Hanning window is applied, and its real part is plotted as a function of time. The result is a discrete short-time 
Fourier transform. A magnitude and time-frequency representation of the sliding DFT is commonly referred to as a 
spectrogram. As shown, the results seem comparable, but the sliding DFT suffers from increased smearing in the 
resolution of the harmonics (lower left comer) and at the end of the sweep. It is well known that the inverse sliding 
transform suffers from loss of resolution, or “time smearing,” 13 and is also computationally inefficient. Furthermore, 
the short-time Fourier transform suffers from fixed frequency resolution for all times; whereas, the wavelet filter 
bank approach allows adjustments in time and frequency resolution pertaining to the signal characteristics. 

APPLICATION TO F/A-18 SRA AEROELASTIC DATA 

A procedure for deriving transfer functions from the SRA wingtip excitation data is described. Time-frequency 
analyses are applied to this data for enhancing transfer functions and improving system identification methods such 
as subspace identification. 


Transfer Function Derivation 

Estimation of transfer functions from the SRA data was predominantly used to identify state-space matrices for 
prediction of flutter boundaries. Because the input was corrupted by exciter vane anomalies and inconsistencies, 
extracting information relevant to linear system identification and considering the rest as being detrimental to linear 
stability estimation was necessary. Therefore, assuming the system to be identified is linear, the extraneous energy 
at frequencies other than the excitation frequency is considered as corruptible unmodeled dynamics, noise, or both. 
Along the main harmonic of the input will also be noise (as with output), but the signal-to-noise ratio is assumed to 
be sufficiently high to mask the majority of the noise corruption. 





(a) Classical Fourier transfer functions. 
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(b) Corresponding time-frequency generated transfer function. 

Figure 5. Transfer functions computed with classical Fourier technique and time-frequency method. Left exciter 
input to accelerometer output (left column). Right exciter input to accelerometer output (right column). 
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Signal in time/frequency domain 


Signal in time/frequency domain 



(a) Morlet continuous wavelet transform. (b) Discrete short time Fourier transform (right) 

Figure 6. Morlet continuous wavelet transform compared to the discrete short-time Fourier transform of SRA 
exciter data. 


The procedure for SRA transfer function derivation from the symmetric and antisymmetric exciter inputs 
is as follows: 

• Compute the Morlet- filtered CWT coefficients of the inputs n,(co,) and « 2 (t o-) from each of the two exciter 
vanes (operating simultaneously for symmetric u s and antisymmetric u° sweeps) at each filtered center 
frequency co ; . 

• Compute y(C0j) for each output y' and y“ at each filtered center frequency to,. 

• Estimate transfer functions G,(co,) and G 2 (co,) from each exciter vane input to the output: 


• and for each to, 


/(CO,) 


m/((0,) «/( CO,.) 

Gjfto,) 

/(CO,) 


m/co,) m/( CD,-) 


GJ (co,) 


«/(C0,.) m/((0,.) 

-1 

/(C Of) 

G 2 (co,) 


«/(©,■) n/(co,) 


/(CO,.) 


(14) 


(15) 


Time-Frequency Analysis 


A demonstration of applications using wavelet time-frequency maps, called scalograms, for data enhancement 
in system identification of SRA structural dynamics follows. In figure 7(a), a scalogram showing the original input 
signal of a linear frequency sweep from the exciter is displayed. A primary harmonic is clearly seen as well as other 
dynamics and noise. The main harmonic is the desired part of the signal because it tracks the programmed input to 
the exciter. The actual input to the exciter is much more complicated due to nonlinear aerodynamics and distributed 
loading across the exciter vane. By creating a mask around the main harmonic CWT of the input and applying it 
to the CWT of the output signal (accelerometer, fig. 8), the dynamic properties of the system can be identified. 
Figure 7(b) represents the masked input CWT. 
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Signal in time/frequency domain 


Signal in time/frequency domain 
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(a) Scalogram of original SR A exciter signal. 

Figure 7. Morlet scalograms of original input and cleaned input. 


(b) Scalogram of cleaned SRA exciter signal by 
feature filtering the CWT. 


Signal In time/frequency domain 
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Signal in time/frequency domain 
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(a) Scalogram of original accelerometer response. 


(b) Scalogram of cleaned accelerometer response. 


Figure 8. Morlet scalograms of raw and cleaned accelerometer output. 


Signal reconstruction uses the redundancy of the wavelets to estimate a time signal which the CWT best 
approximates. The real transform Hj>( co) assures an in-phase reconstruction compared with the original signal. 
Figure 9 shows the original (unprocessed) and reconstructed clean (processed) signals from figure 7. The same 
procedure is performed on left and right exciter inputs and an output accelerometer signal to generate transfer 
functions in figure 10. Note the extreme contrast between the original results (fig. 10(a)) and the cleaned (fig. 10(b)) 
transfer functions, especially in phase. These plots show the improved noise removal in the cleaned data for 
distinguishing modes below 20 Hz, and the potential for dramatic differences in phase between input and output 
responses upon cleaning. Estimation of a state-space system can then be performed using the enhanced data from 
the cleaned CWT and an algorithm, such as subspace identification. 14-17 
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Reconstructed signal In time domain Reconstructed signal in time domain 



Time, sec g70013 Time, sec g70014 

(a) Reconstructed signal of original accelerometer (b) Reconstructed signal of cleaned accelerometer 

response. response. 

Figure 9. Original signal and cleaned signal reconstruction. 


TF ul -> y - 1 TF u2 -> y - 1 Left exciter Input Right exciter Input 
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(a) Transfer functions using raw time-frequency data. 
Left exciter (left column) and right exciter (right 
column) to accelerometer output. 
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(b) Transfer functions using cleaned time-frequency 
data. Left exciter (left column) and right exciter (right 
column) to accelerometer output. 


Figure 10. Transfer functions using raw and cleaned time-frequency data. 


An example of filtering the undesirable features of complicated input-output signals using the time-frequency 
representation will be a double logarithmic sweep from the SRA excitation system. Figures 1 1 and 12 show planar 
and three-dimensional scalograms of the original and desired sweeps. Harmonics from the strain gauge input 
measurement can be readily detected in figures 11(a) and 12(a). The harmonics indicate nonlinear exciter vane 
response from the rotating slotted cylinders at the wingtips, which is deemed undesirable for subsequent linear state- 
space identification methods. Therefore, the input signal is modified by extracting the desired time-frequency map 
from the scalogram (figs. 1 1(b) and 12(b)) and reconstructing the time domain input signal. 

As exemplified, general maps in the time-frequency plane can be effectively filtered. Also, such signals would 
be very difficult to threshold because the harmonic wavelet coefficients are of comparable amplitude. Hence, a 
method of automating this feature filtering process will generally not benefit from typical thresholding schemes. 18 
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Signal in time/frequency domain 



(a) Original input. 


(b) Cleaned input. 


Figure 1 1 . Scalograms of original double-sweep logarithmic input and cleaned input. 


3D representation of time/frequency decomposition 


3D representation of time/frequency decomposition 



(a) Original input. (b) Cleaned input. 

Figure 12. Three-dimensional scalograms of original double-sweep logarithmic input and cleaned input. 


Subspace System Identification 

Estimated transfer functions with a classical Fourier method and the time-frequency localization approach are 
compared in figure 5. The cleaned transfer functions are clearly superior in illuminating the modal responses relative 
to the noise. This distinct difference between the two approaches would considerably influence a system 
identification procedure. 

Example plots in figure 13 compare SRA transfer functions from a typical aeroelastic data set with transfer 
functions generated from the frequency domain algorithm 15 outlined in the appendix. Transfer functions are from 
exciter inputs to an accelerometer output and are derived from unprocessed and processed time-frequency data. 
Matches between unprocessed data and the subspace estimate are in figure 1 3(a), and figure 1 3(b) matches processed 
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(a) Transfer functions derived from raw time-frequency data. 




(b) Transfer functions derived from cleaned time-frequency data. 

Figure 13. Discrete-time frequency domain subspace identification results. Solid lines are transfer functions from 
exciter input to accelerometer output computed with raw (a) and cleaned (b) time-frequency data. Dashed lines are 
transfer function estimates from a frequency domain subspace identification algorithm processing the solid line data. 
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data with subspace estimates. The solid line on these plots are transfer functions derived from the time-frequency 
procedure, and the dashed lines are subspace identification estimates. Agreement is consistently better between 
processed transfer functions and the identification, compared to the original transfer functions, thereby providing an 
increased degree of confidence in the system quadruple estimates when using the feature-filtered data. 

APPLICATION TO F/A-18 HARV AEROSERV OEL ASTIC DATA 

The excitation signals used for generating aeroservoelastic excitation data are described. Morlet decompositions 
of data from these maneuvers are used for analysis. 

Excitation Mechanisms 


Schroeder-phased harmonic signals were compared to sine-pulses \ — — J because linear and logarithmic 

sweeps are deficient at low and high frequencies, respectively. OBES Schroeder-phased and sine-pulsed signals 
were chosen to excite the HARV at the same range of flight conditions of 1 g, 5° to 70° AOA and 30,000 ft altitude. 
Analysis of the Schroeder-phased data compared to linear sweep and sine-pulse data at the same conditions is 
reported in reference 19. 

Figure 14 shows a Schroeder wave input to the aileron and the resulting aileron position and roll rate. Figure 1 5 
shows a sine -pulse input to the aileron with position and roll rate. These signals were input to the surfaces 
sequentially in the following order within 70 sec: symmetric stabilator, differential stabilator, ailerons, rudders, pitch 
thrust vector vanes, and yaw thrust vector vanes. 

Decomposition of Schroeder-Phased Harmonic Data 

Figure 16(a) shows an example Morlet scalogram of a 5 to 25 Hz Schroeder input into the pitch thrust- vectoring 
command at 10° AOA. The nature of the time-frequency response is the sum of two distinct sweeps, one sweeping 
the low frequencies followed by the high frequency sweep with little overlap. All frequencies from 5 to 25 Hz are 
being excited equally. The time-frequency-filtering procedure used to analyze the SRA data is applied in 
figure 16(b), and subsequent analysis applies accordingly. In contrast to the SRA data, the OBES inputs are 
relatively clean because they are generated digitally by the flight system. Figure 17 shows the corresponding transfer 
functions of structural response from the OBES inputs of pitch and yaw thrust vectoring. The procedure enhances 
the modes near 15 Hz and above relative to the other frequencies, which is justified from figure 1 8 where the output 
accelerometer response scalogram is plotted. The response is dominated by energy around 1 5 Hz and above relative 
to other frequencies, with little response below 15 Hz. Note the order of magnitude difference. 

Decomposition of Sine-Pulse Data 

Figure 19 shows scalograms of sine pulses. In contrast to the Schroeder-phased signal, sine functions contain 
the majority of their energy in a small amount of time. It seems that processing of these data is unnecessary because 
nearly all the energy is concentrated in a well-defined area in the time-frequency plane localized in time. Figure 20 
shows corresponding transfer functions of structural response from the OBES inputs of pitch and yaw thrust 
vectoring. Noise in the spectral response is drastically reduced in figure 20(b) compared to figure 20(a). In 
figure 20(a), some difficulty exists in discriminating modal response from spurious noise peaks; in figure 20(b), the 
modal response is more distinct. These transfer functions demonstrate that considerable improvement in data quality 
is achieved by time-frequency filtering. Time localization of the signal degrades the frequency response unless this 
property is exploited in the analysis. In figures 19(b) and 20(b), the time-localized energy is extracted from the 
input-output signals to define a much cleaner frequency response from which to perform system identification. Note 
the order of magnitude difference in gain between raw and cleaned magnitude responses. Also note that the 
responses for the same flight condition at 10° AOA and same input-output signals are very different between the 
two types of excitation, Schroeder-phased and sine-pulsed, either unprocessed or processed. 
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Figure 15. Example 
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(a) Raw scalogram. 


(b) Cleaned scalogram. 


Figure 16. Morlet scalograms of Schroeder-phased harmonic input. 
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(a) Transfer functions derived from raw time- 
frequency data (left). Pitch thrust vector input (left 
column) and yaw thrust vector input (right column). 


(b) Transfer functions derived from cleaned time- 
frequency data (right). Pitch thrust vector input (left 
column) and yaw thrust vector input (right column). 


Figure 17. Transfer functions from raw and cleaned Schroeder-phased data. 
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Figure 18. Accelerometer response sc 
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(a) Raw scalogram. 


Figure 19. Morlet scalo, 
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(b) Cleaned scalogram. 
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(a) Transfer functions derived from raw time- 
frequency data. Pitch thrust vector input (left column) 
and yaw thrust vector input (right column). 
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(b) Transfer functions derived from cleaned time- 
frequency data. Pitch thrust vector input (left column) 
and yaw thrust vector input (right column). 


Figure 20. Transfer functions from raw and cleaned sine-pulsed data. 


CONCLUSIONS 


Multiresolution analyses of flight vibration data on two F/A-18 aircraft were performed with nonorthogonal, 
nondilated Morlet wavelets. This procedure offered enhanced capabilities for time-frequency feature extraction and 
filtering and for subsequent system identification. Such a time-frequency approach allows improved visualization 
and understanding of the signal information content. Noise and other distortion dynamics can be identified, reduced, 
or eliminated. Automated thresholding would be difficult in the majority of situations when distortion dynamics are 
significant relative to the signal content. This proposed approach succeeds when standard threshold techniques are 
inappropriate. Identification schemes used to extract modal data, state-space representations, or stability boundaries 
demonstrate improved performance with these procedures. 
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APPENDIX 


An outline of system identification is presented in time and frequency domains. Proposed advantages and 
experiences with algorithms using structural excitation data is described. 

Subspace Identification of System Quadruple 

Subspace identification methods attempt to find a linear, time-invariant, state-space realization 

x k+\ = Ax k + Bu k + (a k (16) 

Vk = Cx k + Du k + V k (17) 

based on the input vector sequence u k and output sequence y k . The state vector x k is not directly observed, so it is 
only identifiable to within a similarity transformation. White noise vectors 0) k and v k are zero-mean process and 
measurement noise, respectively. The state-space model identification problem is then as follows: 

Given N input-output samples {u k , y k ) of the unknown system, find a consistent estimate of the 
state-space quadruple [ A,B,C,D ] up to a similarity transformation T. Determine the Kalman gain K 
such that the output of the estimated system 

**+1 * Ax k + Bu k + K(y k -Cx k ) (18) 

y k = Cx k + Du k (19) 

is the minimum variance one-step ahead prediction of the recorded output. 

Subspace algorithms are regarded as powerful modem techniques for linear system identification because of 
numerous advantages. For example, 

• They are noniterative, nonparametric, and linear. 

• They are numerically robust algorithms. 

• They produce nonparametric state-space models. 

• Their MIMO and SISO algorithms are identical. 

• They have both time and frequency domain versions. 

• They have model insensitivity to matrix perturbations. 

• There is good performance for identification of high-order (structural) systems. 

For these reasons, the following algorithms are outlined with emphasis on experiences in performance, using data 
processed with Morlet wavelet filtering. 

Block Hankel matrices that play an important role in subspace methods are 


U = 



u | u 2 

u 2 u 3 

u k u k + 1 
u k+ 1 u k + 2 

u 2k u 2k + 1 


' U N - 2k + 1 
' u N-2k + 2 

■ u N-k 
• U N -k + 1 

U N 


( 20 ) 
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where U and Uy designate past and future samples. A similar construction is made for the output matrix Y. The 
extendec/observability matrix is defined as 


T = 


c 

CA 


CA 


s - 


1 


( 21 ) 


Fundamental to subspace algorithms is the following matrix equation 

Y = TX + HU + V (22) 

where Y, U, and V are Hankel matrices formed with the output samples, input samples, and noise, and H is the 
Toeplitz matrix of Markov parameters 



D 

0 

0 

. 0 


CB 

D 

0 

. 0 

H = 

CAB 

CB 

D 

. 0 


CA r 2 

B CA S ~ 3 

c _ 4 

B CA B . 

D 


The following types of subspace identification methods were compared for this study: strictly time domain, 16 
combined time and frequency domain," and frequency domain. 14 ’ 15 - 20 Time domain algorithms had the worst 
performance, especially in extracting B and D matrices, so a different approach was proposed. The A and C matrices 
were derived from the time domain data, then the B and D matrices were computed from matching the state- space 
model with estimated transfer functions. This approach gave better results than using a strictly time domain version 
in getting the B and D matrices, especially when using the filtered transfer functions from the Morlet time-frequency 
method. 11 Finally, two algorithms using frequency domain data were tried to take advantage of the improved transfer 
functions. All algorithms will be outlined here for completeness. 

Time Domain Algorithm 


Because the transfer function estimates are based on at least two distinct data sequences from the SRA data 
(symmetric and antisymmetric), a method was used to determine system A and C matrices from multiple data sets. 21 
Briefly, A and C are derived from the following technique (C is p by n). 

• Find P such that P = TQ and rank(P) = rank(T)- In the noise-free case ( V= 0), P = YU 1 because HUU 1 = 0 
in equation (22) (1 denotes an orthogonal complement). For two data sets, P = [Pj P 2 ] = HQ i Q?\- The 
following steps remain the same. 


Decompose P = USV, where 5 = 


Evaluate: A = u\u { where respectively omit the first (last) p rows of U j, and t denotes 

pseudoinverse. C is the first block of U-y . 

Premultiply (2) by T 1 , so T X T = 0 and post-multiply by U + leading to a least squares solution of 


' Sj 0 ' 

, and U = 

i \ 

V 1 

< 0 0 , 


u 2 

V 1 ) 


are rank(P). 
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Combined Time and Frequency Domain Algorithm 

In the combined time-frequency domain algorithm, the B and D matrices are found from discrete transfer 
function G(zp, where z k = e ^ <£>k and = f° r M data, k = 0 , ... , M 


G(z k ) 



(24) 


for all frequencies co^, the least squares solution separating real and imaginary is 




1 

0 

1 

£ 

W .L 

i 

t 





G( Z] ) 

B 

D 

~ 

C(z 2 /-A) 1 I 


G(z 2 ) 



C(z m I-A )~ 1 I 


°(z m ) 


(25) 


Frequency Domain Algorithm 


Frequency domain algorithms 1415 - 20 are similar in context. In particular, one method 15 starts with the premise 
that the transfer function G of the system has a real-valued impulse response. As a result, data on [0,7t] can be 

extended to [n,2n] by taking the complex conjugate and continuing as follows. 

* 

* Extend transfer function samples to the full unit circle: G M + k := G M k - 1, .. . , Af- 1 . 


1 — ,2Af-l 

Define the real block Hankel matrix: 2A/-point inverse DFT /z,- = — - > G k e 

f— -I 2/Vi k — 0 


jlnik 
2 M 


H : = 


hq + r- 


where {q 9 r) > n and q + r < 2M. 


• Compute the SVD: H = UlV . 

• Determine system order n by inspecting singular values and partitioning: H = [U S U 0 ] 

• Determine system matrices A and C as A - (U s )^ U s and C is the first /7-block of U s . 


Is 0 
0 to\ 


v T 

v o 


• Solve least squares, as above, for B and D . 


Other algorithms use auto- and cross-spectral data in the same framework 14 or identify a Laplace domain 
[A,B y C y D] from discrete transfer function data co^ using numerically well-conditioned Forsythe polynomials. 20 
These frequency domain algorithms perform well for high-order structural systems. 

To recap, even with processed data from the time-frequency filtering procedure, the subspace time domain 
methods perform poorly with rapidly varying inputs. Including some frequency domain information in the subspace 
identification method will improve the estimates. Furthermore, enhancements are necessary in data processing for 
frequency domain data without assumptions of stationariness, time-invariance, periodicity, and white noise. Time- 
frequency filtering methods offer a solution. Other examples using a variety of estimation algorithms can be found 
in reference 22. 
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